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This paper analyzes dynamical properties of small Fermi-Hubbard and Bose-Hubbard systems, 
focusing on the structure of the underlying Hilbert space. We evaluate time-dependent quantities 
such as the return probability to the initial state and the spin imbalance of spin-1/2 fermions. The 
results are compared with recent experimental observations in ultracold gases. For the symmetric 
two-site Fermi-Hubbard model we find that the spin imbalance and the return probability are 
controlled by two and three frequencies, respectively. The spin imbalance and the return probability 
are identical for the asymmetric Falicov-Kimball limit and controlled by only one frequency. In 
general, the transition probabilities between the initial state and energy eigenstates depend strongly 
on the particle-particle interaction. This is discussed for "self trapping" of spinless bosons in a 
double-well potential. We observe that the available Hilbert space is reduced significantly by strong 
interaction. 

PACS numbers: 31.15.ac, 34.50.-s, 03.67.Lx 

I. INTRODUCTION 

The dynamics of many-body quantum states has been studied with high accuracy in ultracold gases, 
either for bosons [H, S 13) 0] or for fermions [1, 0| • The main advantage of employing an ultracold gas 
is that there are many parameters in the experiment that can be adjusted to control the initial state as 
well as the dynamics of the system. This means in particular that an ultracold gas can be prepared in 
almost any state |\E'o), not only in the ground state. After its preparation at time t = 0, the state evolves 
in time t > Q for a system with Hamiltonian H , which describes the kinematics and the interaction of 
the atomic particles, according to 

\^t)=e-'"'\^^) . (1) 

In the following we will consider the simplest possible case, namely a model with two sites. As a first 
example we study two fermions with spin 1/2, a second example is a system of N spinless bosons in a 
double well. Despite of their simplicity, these examples reveal dynamical properties that may have also 
implications to more complex systems. For instance, an interesting question is whether or not all states 
of the underlying Hilbert space can be reached from a given initial state, which time scales appear and 
how this can be controlled by adjusting the parameters of the system such as the tunneling rate and the 
interaction strength. 

The dynamics of atomic systems with a finite number of atoms has been an active field of research 
for several decades, with a focus on oscillating behavior between different atomic states [1, |^. More 
recently, ultracold gases in an optical lattice have been a successful platform for such investigations 
P, 0, [lO]- In particular, small systems of spin 1/2 fermions have also attracted considerable attention 
by the quantum computation community [a. Illl|. A typical initial state is a (paramagnetic) Mott state, 
where the orientation of the individual spins can be used for information storage. Since for N particles 
this state has an exponential degeneracy 2''^, the time evolution in terms of, for instance, a Hubbard 
Hamiltonian can lead to a complex dynamics. Of particular interest is how the spin population of a given 
site changes with time. 

Experiments with a Bose gas in an optical double well have revealed that the population dynamics is 
controlled by direct tunneling of atoms and by a second-order (superexchange) process, characterized by 
two frequencies P, [l3] • The appearance of the characteristic frequencies in finite bosonic systems was 
calculated in strong-coupling perturbation theory [l^, [l^ and in time-dependent Hartree-Fock approxi- 
mation [3, [l^. In a more recent paper by Trotzky et al. [3| the spin imbalance of two spin-1/2 atoms in 
a double well was studied. Such a system can be described by a two-site Bose-Hubbard model, which is 
often called the two- mode approximation [l^ . With increasing values of J /U , where J is the tunneling 
rate and U is the local interaction strength, the oscillating spin imbalance was increasingly influenced by 
a second frequency. This implies that both, the tunneling rate of the atoms as well as the interaction 
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strength, determine the dynamics of the many-body system. In this paper we wiU study this effect in 
terms of a Fermi-Hubbard modeL The focus wiU be on the return probabihty of the many-body state. 
This is an important quantity for recovering information which was stored in the initial state. 

There are several options for an analytic calculation of physical quantities in a finite many-body 
system. Perturbative methods are very successful and reliable approaches to physical problems and can 
be considered as exact if their regimes of validity are respected. On the other hand, they fail in most 
cases near a singularity, where their validity is violated. Nevertheless, they can still be used in these 
cases as approximation methods by using an asymptotic partial summation of the perturbation series. 
Hartree-Fock approximations^ belong to this type of approximations and have been used in the case of 
finite many-body systems [ij, • Other self-consistent approaches to many-body systems have been 
very successful, such as the dynamical mean- field theory [l^ . A disadvantage of all these approaches is 
that the (nonlinear) self-consistent equations are very complex, such that their treatment usually requires 
intensive numerical work. 

An alternative to the perturbative approach and to self-consistent approximations is the recursive 
projection method (RPM) [13, [III- It is a systematic exploration of the Hilbert space, using a recursive 
calculation of the resolvent (z — H)~^. The latter can be obtained from the Laplace transformation of 
Eq. Il]) (cf. Sect. IIII Ap . The RPM enables us to extract systematically the poles of the resolvent in a 
subspace of the underlying physical Hilbert space. This method, combined with a truncation of recursion, 
is related to the Lanczos procedure [l^ . The RPM has been explained elsewhere in the literature [ill > 
here we give only a brief summary and apply it to the Hubbard models in Sects. IIV[ fVl 

The paper is organized as follows: After a short description of the RPM in Sect. |TT] we define the 
Fermi-Hubbard and the Bose-Hubbard model, the many-body return probability, the spin imbalance and 
the many-body spectral density in Sect. IIIII In Sect. IIII Al the dynamics is discussed in terms of the 
resolvent. Then the recursive projection method is applied to the two-site Fermi-Hubbard (Sect. IIV[) and 
to the two-site Bose-Hubbard model (Sect. |Vl to derive dynamical properties of these systems. Finally, 
the results of the RPM are discussed in detail in Sect. IVII In particular, frequencies and transition 
probabilities are calculated and their dependence on the model parameters are discussed. 

II. RECURSIVE PROJECTION METHOD 

The structure of our physical system is completely determined by a Hamiltonian H that acts on Hilbert 
space H. Then the central idea of the RPM is that the dynamics starts from an initial state which lives 
in a subspace Ho C Ti. It should consist of a basis that is dynamically separable, meaning that the 
Hamiltonian does no allow to move directly from one basis state to another. Once this subspace has been 
chosen specifically, the rest of the RPM is entirely determined by the Hamiltonian H, and the dynamics, 
given by the time-evolution operator, decides which part of Ti is relevant. This depends on the energies 
associated with the remaining Hilbert space and on the transition probabilities. Now we project with 
projector Pq onto the Hilbert space Hq. The corresponding projected resolvent is Go{z) = Po{z — H)^^Po. 
Then the RPM includes two steps: 

(1) Hilbert space 7^2^+2 (j — 0, 1, ...) is created by acting with the operator {1 — Pq — P2 — ■ ■ ■ — P2j)HP2j 
on 7i. In other words, a basis set from all the states created by (1 — Pq — ^2 — • • • — P2j)HP2j is a basis 

of 7^2^+2 ■ 

(2) evaluating the resolvent G2j on TL2j through the recurrence relation 

G2, = (z-if^,);; (2) 
with the effective Hamiltonian on 7i2j- 

H2j = P2j HP2J + P2j HG2J+2HP2, . (3) 

remarks: (I) It should be noticed that the construction of the sequence 7i2j (j — 1,2,...) implies that 
Hilbert space TL2j is orthogonal to 7i2j' for j' ^ j. Moreover, the recurrence relation never returns to 
previously visited subspaces. This can be represented schematically as a Russian-doll structure, shown in 
Fig- Q] (II) The recursion terminates in a finite dimensional Hilbert space with the effective Hamiltonian 



ff^„ - P2nHP2. 
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FIG. 1: The schematic structure of the sequence of Hilbert spaces 7^2 j for j — 0, 1, created from the initial 
Hilbert space Hq by applying the Hamiltonian H, resembles a Russian doll. The resolvent Go is projected onto 
Tlo and the outer shells contribute as poles to Go (cf. Sect. |ll|. The recursion removes one shell (mother doll) 
after the other until one reaches the baby doll Ti.o- 



This is the only effective Hamiltonian that is exphcitly given, provided we know the projection P„. In 
order to use this as the initial effective Hamiltonian, we introduce k = n — j as the running index in the 
recurrence relation. Then we have 

G2(n~k) ^(z- H2(n-k)) > (4) 

and with = G2(n-k), ^fc = ^2(n-fc) obtain the recurrence relation 

5fc = (2 - /»fe)2(Lfc) . (5) 

where 

hk = P2(n-k)HP2(n-k) + P2(n-k)H <Jk-lH P2(n-k) , = P2nHP2n ■ (6) 

This means in terms of the Russian doll that, in order to evaluate Go (i.e. the baby doll), we must inherit 
the properties of all generations of mother dolls G2j (j = 1, n), using th recurrence relation iteratively. 



III. MODEL 



The Fermi- Hubbard (FH) model describes locally interacting fermions with spin a =t,i, the Bose- 
Hubbard (BH) model locally interacting spinless bosons. It is defined by the Hamiltonian H = Hj + Hj, 
where Hj is the tunneling term between the sites 1 and 2. For fermions this reads 

= - ^ Ja{c\^C2a + cl^Ci^) + U{n^nii+n2^n2i) (7) 

with fermionic creation operators cj^, annihilation operators cjo-, and density operators nj„ = '^ja'^i'^- 
There are two special cases, the symmetric FH model with J| — J| and the Falicov-Kimball limit J| > 0, 

Jt = [13. 

For bosons with bosonic creation operators bj, annihilation operators 6j, and density operators nj = fej fej 
the Hamiltonian reads 

H = -J{b\^b2a + btbi,) + U[ni + nl) . (8) 

This Hamiltonian is also known as the two-mode approximation of a continuous Bose gas in a double- well 
potential (l3 |. 

For fermions the local (Hubbard) interaction Hj = [/(ni|ni| + ^21^21) can be diagonalized with 
product (Fock) states 

1(71,(72) ((7j =0,T,i,Ti) , (9) 
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where this product state is a paramagnetic Mott state in the restricted case ctj =t,i (singly occupied 
sites). Any eigenvalue of Hi with respect to a Mott state is zero. 

The interaction term of the BH model reads Hj = U{n\ -\- n\), whose eigenstates are also product 
states 712) [rij > 0) with eigenvalues U{n\ + n^). The tunneling term gives for these states 



Hj\ni,n2) = - J{h\bi + b\h2)\ni, 



n2) 



= - J\J {n2 + 1)^1 |ni - 1, ^2 + 1) - J 1/^2(^11 + 1)1^1 + 1, ^2 - 1) • 

States with nj — 1 < do not exist. 

In general, the eigenstates \Ej) of a Hamiltonian H with energy Ej 

H\E,)=E,\E,) 

can be used to describe the time evolution of an initial state |^'o) at time i = to the state |^'t) at later 
time i > by 

= e-»^*|*o> = ^e-»^^*|i?,)(i?,|vI/o) . (10) 
i 

There are two interesting quantities directly related to this expression. One is the return probability to 
the initial state P{t) — |(5'o|^'t)P, the other is the spin imbalance between the two sites in a double well 
i 

Ni,2{t) = ^(*t|"Ti - "ii + ^12 - nT2|^'t) ■ (11) 

The latter describes the exchange dynamics of the two spins f and j, located at the two sites. 7Vi_2 
vanishes if the sites are not singly occupied. 

A. The Resolvent 

A Laplace transformation of for positive time gives the resolvent that acts on the initial state: 

/>oo poo 

-I e'''\^t)dt = -i e''''e-'"'dt\^Q) ^{z~ H)-^\^a) (/mz > 0) . (12) 
JO Jo 

Or equivalently, we can express the time evolution operator in terms of the resolvent [z — H)~^ as 

e-^*= /(z-i/)-ie--*^, (13) 

where the closed contour F encloses all (real) eigenvalues of H. The return probability to the initial state 
is obtained from the inverse Laplace transform of the resolvent through the relation 

(*o|*t)= j^e-''\^o\{z-H)-^\^^)—. (14) 

The spectral representation of Eq. pH]) allows us to write 

^ j 3 

This means that the matrix element of the resolvent 

(*ol(.-i^)-^l^o)=E^f^ (16) 
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is a rational function with poles z = {j ~ 0, 1, ...). The many-body spectral density can be calculated 
from this expression for z = E + ie as 

j \ J ) 

By plotting this expression as a function of E, we can identify graphically the poles Ej [j — 0, 1, ...) of 
Go and the overlap of \Ej) with the initial state. The energy levels Ej are the locations of the Lorentzian 
peaks and \{Ej\^Q)\'^ / e correspond to the height of the Lorentzian peaks. The knowledge of Ej and 
|(£'j|5'o)P enables us to determine the expression for the dynamical overlap function in Eq. ([15]). 



IV. DYNAMICS OF TWO FERMIONS IN A DOUBLE WELL 



Considering only two fermions with opposite spin, the Hamiltonian acts on a four-dimensional Hilbert 
space and can be diagonalized directly with eigenvalues 0, U,U /2± yx/^/4 + 4J^ for the symmetric case 
J| = J| = J . Also the RPM is simple because it terminates already for n = 1, after creating a single 
pair of empty and doubly occupied sites. The effective Hamiltonians then read 

H'^ ^U, -L-P^H^P^ . (18) 

i^o can also be expressed as a Heisenberg spin Hamiltonian [Tsj : 

H'^^ J2 Po [an {sfS} + SfSf,) + an {S?S^, - 1/4)] Pq, (19) 
<jj'> 

with z-dependent spin-spin coupling coefficients 

«n=24±i£, .„.4^ ,20, 

and with the spin- 1/2 operators 

where a-' (j = x,y,z) are Pauli matrices. The spin components read explicitly 

5- = i(4cx+4cT), 5^ = ^(c|cx-c|ct), 5^ = ^(c|cT-4q) . (21) 

In the symmetric limit J| = J| the effective Hamiltonian in Eq. (|19p becomes an isotropic Heisenberg 
Hamiltonian. The singly occupied eigenstates of Hq are linear combinations of | f, |) and | i, T)- If is 
one of the eigenvalues of (z — U)Hq with A = —{J'\ ± we get for the poles of the projected resolvent 
Go 



(l±Vl-4A/C/2)-{^/t,^/^ (22) 



from the RPM. The asymptotic expressions hold for \/U ~ 0. For strong interaction parameter U, only 
one pole is accessible by perturbation theory, such that the appearance of two poles can be understood 
as a simple non-perturbative effect: the Brillouin- Wiener perturbation theory in powers oi X/U (2ll | gives 
only the low-energy pole for Go(2;), namely z — \/U [I4I, and neglects the high-energy pole z ^ U — X/U . 
We will see later that the low-energy pole is indeed negligible for large values of U . Experimentally, 
however, both energies have been observed in a double well potential 
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V. DYNAMICS OF N BOSONS IN A DOUBLE WELL 



A system of N spinless bosons, distributed over two sites, lives in a + 1-dimensional Hilbert space. 
Using the basis |ni,n2), a special case for the initial state is ni = N, n2 ^ 0: j^I'o) = I-/V, 0). Then all 
projected spaces 7i2j are one dimensional and spanned by |iV — j,j). This leads recurrence relation (cf. 
App. \B) 

1 _ 1 

^ z - U[k^ + {N - kf] -J^{N-k + l)kgk-i ' ^° ^ z-UN^ ' 

Then the projected resolvent for the initial state |^'o) — \N,0) reads 

(7V,0|Go|iV,0) = {N,0\{z-H)-^\N,0)=gN • 

The evaluation of gjv from Eq. is a simple task and leads to a rational function, consisting of a 

polynomial of order iV + 1 in the denominator. 



VI. RESULTS AND DISCUSSION 



Using the expressions derived in Sects. IIVI fVl we evaluate the spin imbalance as well as the return 
probability of the two-site FH model for symmetric tunneling and in the Falicov-Kimball limit. In the 
second part the dynamics of a two-site BH model with TV bosons is studied by evaluating the many-body 
return probability and the many-body spectral density. 



A. Two-site Fermi-Hubbard model 



The spin imbalance of Eq. (jlip can be rewritten in terms of the spin operator 

{^t\St~SI\^t) , (24) 

since = (rij — ni)/2. Then the matrix element for the spin imbalance of the two-site FH model with 
respect to the initial state I^Po) = | T, i) becomes after Laplace transformation (cf. Eq. (fT2|) ) 

(t, i\iz~ K{z))-\Sl - 5|)(z' - H'^{z'))-^\ T, i) = 



2{U - z){U - z'){iU^zz' - AUzz'^ + 2zaU ~ Az'^z'U + ^z^z'"^ - 2z'^a + 2z'aU - 2z''^a + - 6^) 
(4z'2C/2 - %z'^U + 4z'4 + Az'aU - Az''^a + a2 - }P){'iz'^U^ - Sz^ + 4z4 + izaU - iz^a + - b^) 

(25) 

with a = 2(J| + J|) and b — 4J| J|. The denominator of this expression is a product of two fourth-order 
polynomials with respect to z and z', which gives four poles for z and for z' , respectively: 

Zi/2 = U/2 ± v/C/2 + 2a - 26/2 , Z3/4 = U/2 ± ^/U^ + 2a + 2b/2 . 

Since 2a ± 26 = 4(J| ± Jj)^, these poles read in terms of the tunneling rates 

2i/2 = t//2±^C/2/4+(Jt- Jx)2, Z3/4 = [//2±^C/2/4+(Jt+ JJ2 . (26) 

The non-interacting limit U = has the poles Zi/2 — i(>^t ~ "^i) -^3/4 ~ ^i^l + i-'^- they are 
linear in the tunneling rates. 

symmetric HF model: For simplicity, we first consider a symmetric Hubbard model with Jj = J| = J. 
Then the matrix element in Eq. (|25p simplifies because of a = 6 = 4J2, and with J = 1/2 we get for the 
spin imbalance the expression 

(T, i |(z - H[,iz)r\St - Sl)iz' - i/^(z'))-'l T, i) 
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_ 2U'^zz' - 2Uzz'^ + Uz- 2z^z'U + 2z^z'^ - z^ + z'U - z'^ 
~ 4{z'U - z'2 + l)z'z{Uz - z2 + 1) ■ 

It should be noticed here that the poles zi — U have been canceled by the factor {U — z){U — z') in the 
numerator of Eq. (|25p . Thus, only the poles Z2, z^, Z4 contribute to the spin imbalance: 



Z2-0, Z3/4 = U/2 ± v/C/V4 + 1 • (28) 

After transforming back to the time-dependent behavior, the dynamics of the spin imbalance is charac- 
terized by only the two frequencies 23,24: 

{•^t\S^ - S'll**) = ci cos(z3t) + C2 cos(z4t) 

The constant term vanishes because the numerator in Eq. is zero for z = z' — 0. Moreover, the 

coefficients are 

2 2 



+ 4([/ - x/C/2 + 4) x/C/2 + 4(C/ + + 4) 

The spin imbalance as a function of time and the corresponding energy levels are plotted in Fig. [2] for 
U = 2 and two different values of J /U . The coefficient of the higher frequency is substantially smaller than 
that of the lower frequency. Moreover, with increasing J /U the lower frequency as well as the amplitude 
of the higher frequency increases. This behavior is very similar to the experimental observation of the 
spin imbalance by Trotzky et al. 

The return probability is calculated from the matrix element 

(T,i|(z-i?^(z))-i|T,i> 2zc/- 2.2 + 1 



2z(zC/- z2 + 1) 

which again has the poles Z2,z^, Z4 of Eq. (|28p like the spin imbalance. This leads to the time-dependent 
behavior 

where Co = 1/2 and 

1 



C*! — -— , — , C2 — — - 



{u + Vu^ + 4)VlF+4: ([/ - Vc72 + 4)v!72 + 4 

Then the time-dependent behavior of the return probability is characterized by three different frequencies: 
Pt^C^+ Cl + Cl + 2C0C1 cos(z3i) -f 2C0C2 cos(z40 + 2C1C2 cos[(z3 - Zi)t] . (29) 
The ratio of the coefficients C2, Ci is 



Cl ^l + 4/C/2-l ^ ' 

with the asymptotic behavior for U ^ 00. Thus for sufficiently large interaction the oscillating dynamics 
of (^ol^t) is dominated by the lower frequency Z4 = U/2 — whereas for weaker interaction 

two frequencies contribute with similar weig ht, namely Z3/4 = U/2 ± y^WJITT. 

Falicov- Kimball limit: In the case of the asymmetric Falicov-Kimball limit the eigenvalues of Eq. (|26p 
are doubly degenerate for all values of J|. This has some consequence for the dynamics. In particular, 
the spin imbalance of Eq. (|25p becomes 

(T, i l(. - H;^i^))-\St - SDiz' - Hl^iz'))-^l T, 1) ^ (,^,_X'aKlV-U + a) 



and the transition matrix element reads 

(T,i|(z-i/^(z))-^|T,i> 



-2z'^ + 2zU + a 
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Thus the expression of spin imbalance is the product of two transition matrix elements. After transforming 
back to time we get 

with Ci = (z3 - C/)/(z3 - Zi), C2 = -(24 - U)/{z3 - Z4), and = U/2 ± ^U'^/4 + Jf. Then the spin 
imbalance is identical to the return probability: 

(vI/,|5f-5||*t) = |(*o|*t)P = P, , 

where the oscillatory behavior is characterized by a single frequency: 

Pt = Cl + C| + 2C1C2 cos (y'c/2 + 4J2t^ . 

B. Two-site Bose-Hubbard model 

According to the Hartee approximation of the double- well potential , the spectral properties change 
qualitatively when the number of bosons exceeds a critical value Nc ~ J /U , where the regime with 
N > Nc is characterized by "self trapping" . This is a regime in which the system stays in its initial 
state for arbitrarily long times. For N < on the other hand, the regime is characterized by an 
oscillating behavior with frequencies related to the tunneling rate. When the number of particles N 
approaches the critical value Nc, the frequency of the oscillations goes down to zero, indicating a real 
critical behavior. However, this might be an artifact of the classical nonlinear equation obtained by the 
Hartree approximation. For a quantum system on a finite-dimensional Hilbert space we expect no genuine 
critical behavior. Nevertheless, a crossover between two different regimes is possible, where in one regime 
the time scale for escaping from the initial states can be very large and the escape is very unlikely. Such 
a behavior will be studied in the following. 

There are two types of quantities that determine the dynamical behavior on a finite-dimensional Hilbert 
space. These are the energy levels Ek and the transition probabilities |(^'o|i?/c)P between the initial state 
and the eigenstates of the Hamiltonian. Both quantities appear explicitly in the many-body spectral 
density of Eq. (|17|) . For the initial state |iV, 0) the many-body spectral density is plotted for several 
parameter values in Figs. [H [5l This clearly indicates that the distribution of transition probabilities 
is broad for weak interaction, referring to a complex oscillating dynamics, and becomes narrower with 
increasing U . For sufficiently large U all transition probabilities are extremely small except for one (cf. 
Fig. 3]). This implies that the system cannot escape from its initial state. This strong- interaction behavior 
can be linked to the semi-classical self trapping of the Hartree approximation. 



VII. CONCLUSIONS 



The dynamics of finite FH and BH models is characterized by a discrete set of energy levels and non-zero 
overlaps of the eigenfunctions with the initial state. By employing the recursive projection method we 
have calculated these quantities in terms of a many-body spectral function for a double well. Moreover, the 
return probability and the spin imbalance are obtained for the spin-1/2 FH model by the same method. 
Depending on the details of the spin-dependent tunneling process, we have found a single-frequency 
dynamics (Falicov-Kimball limit) and a two- frequency dynamics (symmetric FH model). In particular, 
the spin imbalance of the symmetric FH model with tunneling rate J reveals an oscillating behavior 
with the characteristic frequencies U{1 ± a/1 -I- 1&,P /V^) /2, where the amplitude of the higher frequency 
decreases strongly with increasing interaction U . This behavior agrees well with recent experimental 
observations 

For the double well with N spinless bosons the dynamics indicates a crossover from a many-frequency 
dynamics for weak interaction to a stationary behavior at strong interaction. The latter can be related 
to the self-trapping behavior found in the Hartree approximation. 
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FIG. 2: The four energy levels of the symmetric two-site FH model is plotted as a function of the tunneling rate 
J aX U = 2 (left panel). The double degeneracy at J = is lifted by the tunneling rate J. The spin imbalance 
with initial state l^-o) = | T, i> is plotted for [/ = 2 (right panel) with J/U = 0.05 (full curve), J/U = 0.3 
(dashed curve). There are two frequencies contributing to each curve, namely 03/4 = U{1 ± i/lTl6JVf7^)/2, 
corresponding to the lowest and highest energy level on the left panel. 
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FIG. 3: The six energy levels of the BH model of A'^ = 5 bosons, as a function of the tunneling rate J at 
U = Q.l. At J = there are three doubly degenerate eigenvalues. The degeneracies are lifted by the tunneling 

rate J: The higher the energy, the weaker the lifting of the degeneracies. This behavior is very different from the 
non-interacting case i7 = 0, where all energy levels behave linearly with J: E = ± J, ±3J, ±5J. 
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FIG. 4: 



The spectral density of the BH model of A'^ = 40 bosons, U 

2 , 



0.05, e = 0.02. The height of the peaks 
corresponds to \{Ej\'i'o)\'' /e for the initial state |^'o) (cf. text). For J = 0.05 (full curve) and J = 0.2 (dashed 
curve) there is only one dominant energy level, out of 41 energy levels. This is a signature of "self trapping" due 
to strong interaction {J/U = 1). 
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FIG. 5: The spectral density of the BH model of iV = 40 bosons, U = 0.05, e = 0.02, J = 1 (left panel) and J = 2 
(right panel) . This clearly indicates that with increasing J more and more states can be reached dynamically with 
reasonable probability. Moreover, the spectral density is becoming more symmetric. However, even for J/U = 40 
only about half of the 41 energy levels are visible in the plot. 
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APPENDIX A: RECURRENCE RELATION FOR THE TWO-SITE BOSE-HUBBARD MODEL 

For bosonic operators b and and particle- number states |n) we have 

b\n)=^\n-l), b^n) = Vn+ l\n + 1) . 

Then a recursion step of the RPM the tunneling term Hj changes the number of particles in each well 
by one: 

Hj\ni,n2) = -Jlb\bi + &|62)|"-i, "2) = -J(\/ (r7,2 -|- l)ni|ni - 1,712 + 1) + V"-2("i + + 1,^2 - 1)) • 

Therefore, 7^2^ is two dimensional and spanned by {\ni —j, n2 \ni 712 — j)}, provided rij — j > 0. A 
special case is ni = N, n2 = 0. Then the projected spaces are one dimensional and spanned by \N — 
The projected Hamiltonian gives 

P2jJ{b% + b\b2)\N-j,j) = Jblbi\N~ - JVU + mN^m^j + 1) 

= JVU + ^){N-j)\N - U + + 1) 
and the recurrence relation reads 

{N - J, j\G2j\N - J, j) — . ^ — . 

z- {N - j,j\H'\N - 
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The diagonal matrix elements of the effective Hamiltonians (A'' — j,j\H2j\N — are 
{N - j,j\H!,^\N - j,j) = U[{N - jf + f] + J2(j + 1){N - j){N + l|G2,+2|iV - j - 1, j + 1) 

and 

{o,n\h;,j,\o,n) = un' . 

Using the notation gw-j = {N — j,j\G2j\N — we get from the recurrence relation 

1 



gN- 



z - U[iN - jr+f] - J^j + l)iN- j)gN-i-i 



Finally, wo can use the notation k = N — j which implies j = N — k and the recurrence relation for 
A; = 0,'l,...,A^ 

^ 1 1 . . . 

^'^ z-U[k^ + {N-k)^]-J^{N-k + l)kgk-i' ^° z-UN^ ^ ' 

and 

gN = {N,0\Go\N,0) . 



